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ABSTRACT 

A Multidomain Euler Code has been developed to numerically simulate flows of gases of 
different nature around complex configurations with an emphasis on supersonic and 
hypersonic flows. The main choices concerning computational and numerical aspects are 
described. The code has been written in order to allow easy implementation of new boundary 
treatments or further extensions to more complex sets .of equations. Some typical results 
concerning classical shapes, supersonic and hypersonic Hermes shuttle and transverse hot jets 
ire shown. 

INTRODUCTION 

The development of industrial numerical codes to solve the equations of fluid mechanics 
represents an important investment with some uncertain prospective choices to be made. 
These last years quite an important number of attractive methods have been proposed for the 
simulation of perfect fluid Hows around 3D configurations: finite volume or finite element 
methods, using structured or unstructured grids, having a centered or a non -centered numerical 
scheme, etc... The selection of the most promising method seems to be an hazardous and 
difficult choice, a mutter of compromise between different considerations which are generally 
contradictory. Moreover some important elements of choice, such as the available computer 
technology are difficult to be estimated and it is clear that the evaluation of the algorithm 
efficiency is quite different according to the kind of processing: scalar, vector or parallel. 
From an industrial point of view once the desired level of accuracy Is reached by the method, 
the most important quality of the code is robustness which can lead to select algorithms not 
very computationally efficient. An other important point to be considered is the extension of 
the code to a more complex set of equations of the same family. For instance numerical codes 
solving the Euler set of equations for a perfect gas will probably have to be extended to more 
complex state equations or to multi-species gas sooner or later. The easy implementation of 
new boundary treatments by means of a modular coding, as well as ergonomic considerations 
are also very important matters for future development and use of the code. The aim of this 
paper is tc present a 3D Euler code developed for the numerical simulation of flows of gases 
of different nature keeping in mind all the points stated before. The different gases considered 
are perfect gas, real gas at equilibrium and non-rtactive two-species mixture. In a first pan, 
the choices leading to the architecture of the code are described. The second part deals with the 
numerical scheme and its implementation. Lastly the third part exhibits some first results 
obtained. 


RSF N* 23/1123 AY 


Pag^ 17 


1. ARCHITECTURE OF THE CODE: GENERAL CONSIDERATIONS 

The architecture of the code has been dictated by constraints concerning geometrical 
considerations, computational aspects and the specific nature of the flow. 

ggffli figsal flraidcaiigu 

The treatment of complex geometries has led us to adopt a multi block grid made of several 
structured, possibly overlapping or patched domains. This choice considerably simplifies the 
mesh construction and allows the same generality as unstructured grids. Such multiblock grid 
strategies art currently being developed atONERA [IJ and AEROSPATIALE [2] and will be 
implemented with the code. Another interesting possibility has been introduced to enable 
different kinds of boundary conditions on a given domain face. 

CompunikmAl aspects 

The consiraincs concerning the computational aspects are clearly linked to the available 
computer technology. From this point of view, vector and parallel processing have been 
considered. A natural idea in multidomain codes is to distribute each domain on a processor, 
but this could be not so interesting practically since the domains have very different numbers 
of points in usual situations and most of the processors will have to wait for those dealing 
with the greatest domains. From this synchronization waiting-time point of view a code 
based on the computation of separated planes seems more interesting. Another favorable 
computational aspect of a plane structure is that the working arrays for the numerical scheme 
are addressed by two indices and hence not very expensive in core for the present time 
computers. This cun be very interesting especially in the case of implicit algorithms. Of 
course this way of limiting the arrays indexed by three indices restricts the choice of 3D 
algorithms. For instance ADI m the 3 directions is no more possible. With this choice and 
for meshes of a current industrial size that is to say about 200000 nodes it is possible to 
work without any massive I/O on a computer with a few megi words of memory core. 

Boat saim 

For supersonic (lows it is very important to take advantage of the hyperbolic property of the 
steady Euler equations (see [31 for instance) in the main direction by using space marching 
techniques whenever possible. Doing so the CPU time required decreases by an order of 
magnitude. The fact that a plane of a domain can be computed separately allows to make 
multiblock space -marching computations plane by plane without any additional effort 
Another consequence of this supersonic nature is that for steady problems the nature of the 
computation can be different depending on the domain. A block with a subsonic pocket such 
as the blunt body pun will have to be calculated in an unsteady way, the downstream pan of 
the flow with a space-marching strategy. Possibly some real gas effects will have to be 
considered in some puns of the flow and not in others. The order in which the different blocks 
have to be computed is not obvious, and it is no more possible to iterate similarly m all the 
domains as it can be done for transonic computations. The computation has to be managed 
by a command interpreter to allow a certain flexibility. 

Code orcanimioo 

It results in a code organisation built around 3 key units: a command interpreter which 
assumes the user interface, a plane monitoring unit which decides of the type of the 
compulation, and a plane processor including the numerical scheme. The plane processor is 
described thoroughly in the second pan. We now detail a little more the two first units. 

Command interpreter 

The command interpreter is a language, this means that the input file is interpreted 
dynamically. As a language it has to own classical control instructions such as DO loops or 
IF statements. 


Its main functions are the followings. 
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Monitoring the core : The memory allocation is made so that the amys indexed by the 
whole 3 D data are reduced to the minimum, namely for unsteady applications the three 
coordinates plus the necessary conservative variables (in number equal to that 0/ the 
equations of the system). This kind of data is stored sequentially domain after domain in the 
core by the mean of pointers. Monitoring the core consists in attributing pointers when a 
new domain is created or reorganizing it when a domain has been computed and is no more 
useful. The same kind of monitoring b made for the data necessary to coupling boundary 
conditions where 4 Bilinear interpolating data have » be kept for each node involved [ 4 ]. 

Defininf therwxxtynamic stotts : According to the a erody na mic conditions different kinds of 
systems of units can be appropriate. For low supersonic computations it might be interesting 
to define the variables relatively 10 critical conditions that is to say values defined for a sonic 
flow. For hypersonic flow the use of Mol Her tables leads to S! units. From an ergonomic 
point of view u has also seemed interesting 10 give the opportunity to the user 10 use a large 
set of possible ways of defining thermodynamic states. 

Defining scheme parameters and boundary conditions : For instance ;he Couram number or 
the nature of the boundary condition. 

Initialising domain variables : Many options can be used. Initialization can be made by 
means of formatted or unformatted files, with uniform states previously defined. If necessary 
it can be done on a pan of a domain only. 

Extr act in 1 data This includes not only resuh files to be interpreted but also intermediate 

printings and creauon of PH 1 GS metafiles containing graphical data. 

Collin 1 the plane monitor : And so doing calling the numerical scheme. 

Plane monitoring 

As the calculation is based on the computation of separated planes, different kinds of 
calculations can be made according 10 the way these plana are computed : 

Unsteady flows : The time step is the same for every cell of every domain. 

Pseudo unsteady flows : The solution is advanced in an unsteady way but without taking 
care of the physical meaning of the flow before convergence. 

Space mar chi ns computation : The solution on a physical plane is advanced ull convergence 
before computing the next physical plane. 

Since there are basically two kinds of calculation, space marching and unsteady 
computations, there are also two kinds of multidomain strategies. A plane multidomain 
strategy where all the planes of different domains corresponding 10 a physical plane are 
calculated together, and a more classical way where the solution is advanced on a whole 
domain before considering another one. 


2 . THE NUMERICAL SCHEME 
Formulation 

The unsteady 3 D Euler equations are written In conservation form: 

w t*F x *G y ♦Hj*© 

where, for instance, for a one species perfect gas 
W-Kp, pu. pv, pw.e) 
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F-^pu, p*pu-, puv, puw, (e+p)u) 

G*Xp v . puv. p+pv 2 . pvw. (e*p)v) 

H*Hp w « puw, pvw, p+pw 2 , (e+p)w) 

With pH> 1 Xc- Ulp .(u^v 2 * w 2 )) 

To solve them we use an implicit upwind TVD finite volume scheme of Van Leer MUSCL 
type. To obuin the maximum benefit from the structured organisation of the grid and the 
organisation by plane of the code, the implicit part consists in an ADI like inversion in each 
plane coupled with a G our, -Seidel like relaxation in the third direction. Basically the scheme 
comprises the 3 following steps: 

1) Introduction of a linear distribution for each direction in each cell to compute the cell 
Interfaces 

^i+a. j+X, k*n * jjt * 0 J*ijk * ^ l^ijk * ^ I^ijk 

where U-P^OV) is a set of variables. To preserve stability near discontinuities it is necessary 
to introduce limiters in each direction: 

if aj+i/2- 

g 1 i j k -Iimiier(Q|^i/2.0i.i/2) 

Many limiters have been implemented among them the "minmod" [5],Van Leer [6].Van Albada [7], 
and "supertee" (5) formulations. The set of variables can be chosen among the conservative, the 
primitive (p.u,v,w,p) or the characteristic one. Then the scheme fulfills a monodimcnsional TVD 
property. 

2) Computation of the explicit pan 

AW M p . -<ii/Vol(F^, / 2.F i . 1/ 2-^3j + | / 2-<}j.i/2-Hi c .|/2-H k .i/2> 

where Fj+j/2 (**P- G. H) is an evaluation of the fluxes at an imerfxe of the cell control 
volume by means of an approximate Riemann Solver between the two states on each side of 
the interface calculated in the first step. 

Many approximate Riemann Solvers have been tested: 

For perfect gas: The Van Leer (8], Roe (9J and Osher ( 10] formulations 

For a mixture of non reactive two species gas: Abgrall [11] extension of Roe fluxes and 

Abgrall-Monugnd 1 121 extension of Osher fluxes 

For real gas with an equilibrium assumption: Vmokur-Momagnd [13] extension of Van Leer 
and Roe scheme, and Abgrail-Momagnd [12] extension of Osher scheme. 

3) Computation of the implicit part in each plane 

A^W«AW cxp 

where A can be seen as an approximation of 
(I ♦ At 3F/dw) 

where F is 2 i\ approximation of F x +Gy*H z 

Two approximations have been tested: the linearized conservative implicit formulation of 
Sieger-Warming and the linearized non conservative formulation of Hartcn-Yee [14] or 
Chakravanhy [151. 
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For boundary conditions many treatments have been considered, from Viviand-Veuillot [16] 
compatibility relations to more classical flux treatments and their implicitarions. 

As it can be seen, Van Leer's MUSCL approach presents several advantages : 

The scheme is TVD and so well suited for flows with strong discontinuities. 

The extension to more complicated state equations or multispecies convective flows is 
straightforward by implementing the corresponding new approximate Riemann solver and 
new jacobian of the fluxes. The modifications are made very locally in the computer code and 
the overall architecture is fully kept Moreover all the modifications can be added so that we 
have a unique source for all the different kinds of flows. 

For most of the variants of the scheme, no parameter have to be fixed by the user. 

Programming, .noics 

To implement this numerical scheme the following steps are coded in the plane processor 
which advances the solution on a plane k cf a domain. 

1. Research of intersecting boundary conditions. 

This step determines among the boundary conditions that have been declared those which are 
concerned with this plane. 

2. Treament of geometric singularities. 

In those routines boundary nodes including fictitious points might be moved in order to treat 
some geometrical singularities like axes or some special boundary treatments like those 
dealing with half control volume cell (for instance wall like conditions) or symmetry. 

3. Computation of metrics . 

From the data of either the cell centers or the cell vertex, the attributes of the control volume 
such as volume, diameter of the included sphere or outward normal vectors are computed. 

4. Computation of time step. 

And then repeat the following steps 5 to 11 for the 3 directions i j.k. 

5. Treatment of boundary nodes. 

According to the boundary treatment involved, conservative values of frontier nodes (possibly 
fictitious nodes) might have io be assigned before the computation of the slopes. It concerns 
nearly all the boundary conditions since the fictitious points have at least to be extrapolated. 
It is crucial for matching boundary conditions where the nodes might be interpolated into 
another domain. 

6. Computation of slopes. 

This is the pan 1 ) of the numerical scheme described above, 

7. Treatment of boundary interfaces. 

It concerns boundary treatments for which the values are treated on the interfaces rather than 
at the nodes. 

8. Computation of fluxc 

This is the pun 2) of the numerical scheme. 

9. Treatment of boundary fluxes. 

Here, boundary fluxes are possibly modified. For instance for a wall treatment the flux 
computed by the scheme is replaced by a pressure (lux. 

10. Computation of implicit coefficients. 

This is pan 3 of the numerical scheme. 

1 1. Treatment of boundary implicit coefficients 

According to the boundary conditions the implicit matrices are modified. 

12. Explicit result 

Obtained by Summing the fluxes it constitutes the right hand side of the implicit pan. 

13. Resolution of the implicit system. 

It consists of an ADI tike inversion but several other options such as a Causs-Seidel, Jacobi 
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algorithms can rosily be implemented. 

14, Compatibility relation treatment of a wall boundary. 

It consists in a modification of the values at the frontiers according to compatibility 
relations. 

This decomposition of the algorithm coding allows the implementation of a large variety of 
treatments. To be more concrete, let us take the case of the implementation of an explicit 
wail treatment condition and a matching condition and sec some of their pouiwe 
implementations. 


For a first order like cell minor condition the fictitious node can be imposed symmetrically 
10 the boundary node and so step 5 can be used. . . 

For a second order like cell mirror condition the outward interface value can be imposed 
symetrically to the inner interface value so that step 7 is used. 

For a classical flux treatment, first a half control volume can be defined next to the wall and 
then the wall flux is explicitly computed as a pressure flux. Step 2 and 9 are used in this 

For a compatibility relation treatment of a slip boundary (*ee (161 ) o^y U u usa *- 


For a matching condition among many possibilities the two following ones have been 

implemented; ,. . . ... _. 

for an area with a smooth now. a second order matching condition can be applied. The 
boundary node and its fictitious counterpart are interpolated in the coupling domain. Then 

only step 5 is used. . . . 

for an area with strong gradients a Tint order maiching condition can be used. A constant 
distribution is introduced in the boundary cell for the outward direction and the e *J 4 f ntl 
interface value is imerpolated in the coupling domain. This can be done with steps 3 


3. CODE VALIDATION 

In this pan some preliminary results showing some possibilities of the code are presented. 
The first case permits the validation of the multidomain space-marching strategy. The second 
case exhibits some results obtained on a more complex shape, the Hermes shuttle •" 
supersonic and hypersonic aerodynamic configurations. The last case illustrates a mulubloc 
computation of a non reactive two species gas flow with an adequate refinement. 

Ogive -cylinder-flare configuration 

Experimental results on this configuration are available in (17], In the selected test case the 
Mach number is 2.96 and the incidence 4 degrees. Three space marching computations have 
been made. The first one with a monodomain grid, the second one with a continuous two 
domain grid, the last one with a discontinuous two domain grid (figure 1). Figure » shows 
the Mach number contours on the body and in the plane of symmetry. It can be seen that 
results are very similar except a small wave issued from the intersection of the ogive attached 
shock with the coupling frontier. Figure 3 puts into evidence the very weak influence of the 
different grids on body pressures and shows a rather good comparison of the computed results 
with experimental ones. The small differences between the curves come from the mtersecuon 
of the previous wave with the body. 

Hermes shuulc configuration . . fn 

Two supersonic computations have been nude for a MacH number of and incidences o u 
and 10 decrees. The Mach number contours on the body And in the symmetry plane Art 
displayed on figure 4, comparison of the lift with experimental dau has been found 10 be in 
very good agreement Figure 5 shows density contours obtained on Hermes forcbody for an 
hypersonic flow of Mxh number 12 without Incidence. The fluxes used in this case are the 
Van Leer-Vinokur*Momagne ones. 
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Hot transverse jet calculation 

This simulation deals with the interactions occuring between a jet emerging straight up from 
a flat plate into a flow parallel to this one (18] (see figure 6 for the initial aerodynamic 
conditions). In a practical situation the hot jet has a specific heat rauo different from the 
outflow. This can be simulated numerically by considering non reactive mixtures of two 
species flows. An important experimental feature of these flows is the existence of a 
subsonic area that could not be observed on the initial coarse grid. So in order to capture this 
subsonic pocket a refined domain has been added on in supposed location (see the frame 
figure 7 and the computational grids figure 8). This domain has been initialized with the 
solution obtained in t) ' coarse grid and then the calculation has been achieved only on the 
refined domain (weak coupling). Figure 9 shows the concentration distribution of the hot jet 
in the refined area for the two meshes. One can see the strong influence of the grid refinement 
on the jet shape. Figure 10 shows the Mach number distribution in use flowficld. Clearly the 
subsonic pocket is nicely captured on the refined grid (minimum Mach number of 0.67). 

CONCLUSION 

A 3D multidomain Euler code has been developed, its very modular coding allows the 
implementation of various numerical variants of Van Leer MUSCL scheme and a large 
variety of boundary conditions. The extension to different state equations or multispecies 
flows has been shown to be straightforward. Its ability to handle unpatchcd grids eases 
considerably the multiblock gridding and permits to make refinements wherever wanted. The 
code is monitored by a command interpreter which possesses the necessary flexibility to 
handle supersonic and hypersonic computations around complex shapes. 
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a) Monodomain grid. 



b) Continuous two-domain grid. 



c) Discontinuous two-domain grid. 


Figure 1. Computation grids of the ogive -cylindcr-flare configuration. 
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i) Mooodomain computation. 



b) Two-domain computation 
using a continuous grid 



c) Two^Jomain computation 
using a discontinuous grid. 


Figure 2. Mach number contours on the body and in the plane of symmetry 
for the ogive-cylinder-flare configuration. 
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Figure 3. Computations-experiment comparison. 

Evolution of the cc *.fficient of pressure along the winward side. 
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a) NU.-2J-O-0*- Roe solver b) M.-2.5 - o-10° - Van Leer solver 


Figure 4. Supersonic HERMES computations - Mach number contours on the 
body and in the plane of symmetry. 



Mo*»12 • a»0° - Von Leer- Vinokur- Monugnd solver. 


Figure 5. Hypersonic HERMES forebody computations • Density contours on 
the body. 
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General considerations 

This poster presents a 3D Euler code developed for the numerical simulation 
of flows of gases of different natures around complex configurations with ™ 
emphasis on supersonic and hypersonic flows. The numerical scheme of the MUbLL 
type uses approximate Riemann solvers of the Van Leer, Roe. and Osher types which 
have been developed for perfect gas flows and recently extended io non nzcwt 
mixtures ot two species and real gas flows by Abgrall. Montagne and Vinokur. The 
architecture of the code has been dictated by constraints concerning geometrical 
considerations, computational aspects, the specific nature of the flow, and 
ergonomy. 


Geometrical considerations , ... , . , 

The treatment of complex geometries has led us to adopt a multiblock gnd 
made of several structured, possibly overlapping or patched, domains. This choice 
considerably simplifies the mesh construction. Another interesting possibility has 
been introduced to enable different kinds of boundary conditions on a given domain 
surface. 


Computational aspects . . . . . , — . 

The numerical scheme is based on the computation of isolated planes. I nis 
structure presents several advantages. Working arrays for the numerical scheme are 
plane addressed and hence not very expensive in core for present time computers. 
This can be very interesting especially in case of implicit algorithms. For future 
adaptation to parallel machines, synchronization waiting times are considerably 
decreased since every plane requires about the same amount of CPU, which is not the 
case for domains in practical cases. 


Flow pattern 

For supersonic flows it is very important to take advantage of the hyperbolic 
property of the steady Euler equations in the main direction by using space marching 
techniques whenever possible. Doing so the CPU time decreases by an order of 
macnituded ). The fact that a plane of a domain can be computed separately allows to 
make multiblock space marching computations physical plane by A“* 

without any additional effort The presence of strong shocks in this kind of flow 
leads us to select upwind TVD schemes which have proved to be well suited for 
accurate capture of discontinuities. 
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Ergonomy 

For this kind of application the order in which the domains ire computed is 
not always straightforward. Some domains hive to be converged before others, and 
the type of computation in each domain can be different (pseudo unsteady, space 
marching, time accurate...). Thus, in order to add flexibility and user-friendliness a 
command interpreter has been implemented. 

It results in a code organisation built around 3 key units: a command 
interpreter which assumes the user interface, a plane monitoring unit which decides 
of the type of the computation, and a plane processor including the numerical 
scheme. 

The numerical method 

The unsteady 3D Euler equations are written in conservation law: 

Wj+Fx+Gy+HjxO 

where, for instance, for a one specie perfect gas : 

W»t(ppu.pv.pw,e) 

F* l (pu.p+pu 2 .puv,puw,(e+p)u) 

G* t (pv.puv,p+pv2,pvw,(e+p)v) 

H» l (pw.puw,pvw,p+pw2,(e+p)w) 

with p*(Y-D(e-l/ 2 p(u 2 +v^+w 2 )) 

To solve them we use an implicit upwind TVD-finite-volume scheme of Van 
Leer MUSCL type. To use the maximum potential of the structured organisation of 
the grid and the plane organisation of the code, the implicit part is consumed by an 
ADI like inversion in each plane coupled with a Gauss-Seidel like relaxation in the 
third direction. Basically the scheme is constirued by the 3 following steps: 

1) Introduction of a linear distribution for each direction in each cell to 
compute the cell interfaces: 

Ui+t j+XJc+p-Uij.jt + tgiijit + XgJjjk ♦ pg^ijic 

where U«P*1(W) is a set of variables, to preserve stability near 
discontinuities it is necessary to introduce limiters in each direction, if 


Oi+1/2* Uj+ijjc-Ujjjc 

g i ijk*limiter(ai + l/ 2 .ai-l/ 2 ) 

Many limiters have been implemented among them the "Minmod"(2),"Van 
LeerC3). M Van Albada"(4), and • , Superbee"(2). The set of variables can be chosen 
among the conservative, the primitive (p.u.v.w.p) or the characteristic one. Then the 
scheme assumes monodimensional TVD property. 

2) Computation of the explicit pan 

AWexp- - At / Vo i(Fj+i/2*Fi.i/2 ♦ Gj+j/2 -Gj-l/2+Hk+l/2*Hk-l/2) 

where Fj+j/2 is an evaluation of the fluxes at an interface of the cell control volume. 
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I his is made with ihd help of an 'approximate Riemann Solver between the two states 

on each side of the interface calculated in the first step. 

Many approximate Riemann Solvers have been tested: 

-For perfect gas: The Van Leer(5), Roe(6) and Osheitf) formulations 

-For a mixture of non reactive two species gas: Abgrall(8) extension of Roe 
fluxes and Abgrall-Moniagnd(9) extension of Osher fluxes ' 

•For real gas with an equilibrium assumption: Vmokur-Montagne(IG) 
extension of Van Leer and Roe scheme, and Abgrall-Monugnd(9) extension of 
Osher scheme. 

3) Computation of the implicit part in each plane 
A.AW»AW C xp 

where A can be seen as an approximation of 

(l+Atd£/ dw ) 

where £ is an approximation of r x +Gy+H z 

Two approximations have been tested, the linearized conservative 
implicitation of Steger-Warming or the. linearized non conservative formulation of 
Harten-Yee(ll) or Chakravarthy02). The inversion is made by a DDADI algorithm 
(Diagonal Dominant Alternate Direction Implicit). 

For Boundary Conditions many treatments have been. considered, from 
Viviand-Vcuillot(13) compatibility relations to more classical fluxes treatments and 
their implicitation(14). 

After a description of the organization of the code, some typical test cases are 
presented. First, some results obtained by the implicit multidomain 3D algorithm on 
a classical test case consisting of a- cylinder in a supersonic flow. Then some 3D 
results obtained with the explicit algorithm alone. The interaction of a transverse jet 
with a supersonic flow which enables to see differences obtained when using hot or 
cold jets (14) (y* 1.2 or y»1.4) or when using different fluxes. The spreading of the 
iso-concentration lines which should be gathered on a thin surface theoretically, 
shows the diffusive aspect of the schemes. A computation to demonst/ate the 
possibility to refine locally the mesh by the multidomain capability is shown. Finally 
some flowfields around Hermes at supersonic speeds or its forebody alone at 
hypersonic regime show the ability of the method of dealing with realistic 

geometries. 
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SUMMARY 

A variational method for the optimization and adaptation of structured (rids is prtsented. This 
method is applicable to tiro- or three-dimensional structured mesh as well as to surface mesh. It relies 
upon the introduction of a proper measure of a cell deformation that is derived from basic principles of 
continuum mechanics. Several properties are prescribed, which guarantee the well- posedn ess of the mesh 
optimisation problem and the efficiency of the solution algorithm. The ability of the method for the 
mesh adaptation is also described. Examples of two- and three-dimensional grids are shown and 
illustrate the success of the presented method. 


1. INTRODUCTION 

In recent years, remarkable progress has been observed in Computational Fluid Mechanics; first, they 
have concerned the complexity of the physical phenomena and enable us to consider the solution of 
turbulent flow equations including chemical effects. They have also concerned the improvement of 
numerical algorithms for the computation of approximate solutions so as to obtain more efficiently more 
accurate solutions. Finally, a new field of research • generation and optimization of grids necessary to the 
computations - is now being thoroughly investigated, and even becomes a priority for applied CFD. 

Numerous strategies have been developed and more and more sophisticated algorithms have betn 
introduced for the optimization of mesh, e.g. (l, 2j. In particular, adaptation has recently become a very 
important topic of research. This feature allows for instance a redistribution of a fixed number of nodes 
so as to concentrate the points in regions where large variations of the solution occur, and thus increase 
the accurracy of the solution with a fixed number of nodes. However this localized refinement cannot be 
accomplished in any arbitrary way and the quality of the mesh (regularity, non skewness) must be 
carefully controlled. A method able to meet these requirements would generate the "optimal mesh" that 
would enable us to obtain the best possible solution for a fixed computational effort [3] . In this paper, we 
briefly review a method developed in an attempt to more sharply resolve these questions (see also [4-6]) 
and we show the latest results it has provided. The method realizes a tradeoff between the geometrical 
quality of the mesh, the ability to handle adaptation with respect to a given criterion, and the robustness 
and efficiency of the computational algorithm. 


2. FORMULATION FOR A GRID GENERATOR 
A measure of the grid deformation 

The objective of the method is first to obtain a structured and regular mesh as orthogonal as possible 
in a curvilinearly quadrilateral domain in R 2 or on a surface in R\ or in a hexahedral domain in R 3 ; the 
ability of the method to adapt a mesh is obtained in a simple way as it will be shown. In order to reach 
the objective, we consider the mesh generation as the finite element discretization of a continuous problem 
consisting in finding a transformation x(£) from the reference unit cube (space () into the domain to be 
meshed (space x), and we want to exhibit a proper measure a of this transformation. Four axioms and 
geometrical properties can be stated [4]; they are basic principles of continuum mechanics [7], and they 
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lilow the definition of functional* that are not defined in term* of regularity, orthogonality or cell 
skewness, but rather in terms of measure of the deformation of the current cell (with respect to a 
reference cell for instance); these functionals take into account the whole metric of the grid and thus 
define new criteria that measure the mesh quality. Indeed, it can be shown that functionals satisfying 
these axioms must only depend on the invariants l v l v / 3 of the left Cauchy-Green tensor C of the 
transformation x((): 

<r-<r(/i, /,,/,) (!) 

where 

Jj - tr C ; J 2 » tr Cof C ; / 3 • det C (2) 

and 

F - Vx ; C - F'. F (3) 

Furthermore, we want the orientation of the cell to be taken into account; for this, it is necessary to 
make 9 sensitive to this orientation, that is to say to the sign of the transformation Jacobian /, and 
consider functions depending not only on l v but on /: 

«**(/„/,,/) (<) 

with 

J - det F ; /, - / (5) 

Characterisation of the function of the Invariants 

The next step b to ensure the well- posedn ess of the minimization of these functionals: in order to do 
so, we impose new properties relative to the derivatives of 9 for the so-called rigid body modes which 
preserve the shape and the orientation of the cell. We impose a normalization (c vanish« for any rigid 
body transformation), an equilibrium condition (the measure is stationary for any rigid body 
transformation), and most importantly a convexity condition (the measure is convex around any rigid 
body transformation). These conditions, and in particular the convexity, are the proper mathematical 
properties that ensures the well-posedness of the minimization of 9 as well as the efficient convergence of 
numerical algorithms towards a unique minimum. They also restrict the possible choice for the function 
of the invariants [4j; among the functions satisfying the geometrical axioms and these properties, several 
polynomials can be exhibited, in particular 

'-<?>( /, -2 )+<?,(/, - 2 /,-!)+ < 7 , ( 7 - 1 )* ( 6 ) 

where the constants must verify: 

3 C 3 > 4( C { + C 7 ) > 0 (7) 

Furthermore, it b interesting to point out two classes of polynomials that lead to geometrical 
interpretations. 


• In two dimensions, one considers the function 

■ *U - c ( ! 1 - h - 2 ) + K ( J ~ 1 )* with K > C > 0 

.2 


Thb can be re-written as: 


- C ( /, - 2 J ) + ( K - C ) ( J - 1 f 


( 8 ) 

(9) 


The first term ( /, — 2 J ) can be interpreted when one considers the matrix F and its invariants; we 
have: 


F - 


r*e x » 


y< jr, 


, 3 3 3 3 . 

• A ■ x t + z n +>< + ; J - *( - *, 


( 10 ) 


and therefore: 

A - 2 J - ( z ( - y „) 2 + ( *, + i( ) s (n) 

This term represents a least square formulation of the Cauchy-Riemann relations that ensures the 
conformity of the mesh: 

*( - v n - 0 ; y< + *, - o 


(12) 
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and that can also be written: 


x ( »x,Xk _ (13) 

Furthermore this functionr can be ujed to construct a "quaji-conformal" grid: up to this point, the 
quantity 9 was measuring the deformation of a cell with respect to a square reference cell. This can 
however be generalized and it is possible to consider a rectangular reference cell; this allow* tbe 
prescription of a ratio m h ft" for each cell; generalised Cauchy-Riemaan equations then writes: 

A * *< - * “ 0 1 A * *< + * r e " 0 , . ^ U ) . 

when fi \s k conformity parameter tb*fc netdj to b« introduced in order to htee erutenct of tbe quui- 

cooformsl solution *ad that is computed by: 


I(«/0 ( *! + r! ) d( 

- -2 ( 15 ) 

/ Q (»/«) (*? + *; 

This method allow* the construction of grids with arbitrary refinements on the sides such as the one 
presented on Fig.l which Arina also obtained [8] by directly solving the equations (14, 15). 



Fig.l: Quasi-Conform Mesh with Arbitrary Refinements on the Boundaries. 


The second term in (8), ( / - 1 )*, is interpreted as a least square formulation of a constraint / - 1; 
indeed this penalty term is used as a volume control term that prevents the cell volume V from getting 
too far from a reference value V nf ( / - V / V Ktf ), in particular this term avoids the cell to overlap. 

In two dimensions but on a surface in R*, the functional are still valid, and in particular, the term 
/ — 2 / represents a least square formulation of generalized Cauchy-Riemann relations on surfaces [6j: 

z ( »x,Xn (18) 

where n is the unit vector normal to the surface and defining its orientation. 


In three dimensions, one considers the case where: 
(7j — ■ C and C 5 — 3 C +K 


(IT) 
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3. NUMERICAL EXAMPLES 

The solution algorithm and several resulu that this method has produced have been presented by the 
authors in previous publications. These results showed in particular iu robustness (ability to untangle 
ramdom initializations [4]), its ability to regularize a mesh in the neighborhood of singularities [<J or to 
adapt a mesh with respect to a criterion defined by several physical variables (5, 6]. Here we present two 
examples that illustrate the behavior of the method for the adaptation of 2- and 3-dimensional grids. 

The first example presented show* the ability of the method to adapt a gnd around a physical 
phenomenon. A bow shock develops in the solution of the Euler equations for a supersonic flow past a 
cylinder. A first solution can be computed on an initial mesh (Fig. 3a), and a weight w obtained as i 
function of the gradient of the pressure; the minimization of the constructed adaptivity functional leads 
to a mesh that effectively refines in the bow shock area (Fig. 3b). The solution of the equations on this 
new mesh leads to a much more accurate solution than the one obtained on the initial grid [fi|. 




Fig. 3: Adaptation around a Bow Shock 

A second realistic example is proposed: it regards the three-dimensional adaptation to a bow-$ho<k 
developing in front of a re-entry vehicle in a flow at Mach 5 [9l. The Fig. 4c shows the solution 
computed on an initial 31x20x23 point mesh (Fig. 4a), the mesh is adapted with respect to this solution 
(Fig. 4b). The quality of the solution computed on this adapted mesh has been very much improved 
(Fig. 4d) in particular both the bow-shock and the cockpit- attached shock art more precisely described 
These examples, and more generally the gratifying results that this method have produced, have 
shown iu robustness, iu efficiency and iu ability to improve the solution of complicated flow problems 
we think that these features are due to a sound mathematical basis which is necessary to achieve progress 
in the development of grid generation techniques. 
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